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ABSTRACT 

Context. It has been suggested that convection in Red Supergiant (RSG) stars gives rise to large-scale granules causing observable 
surface inhomogeneities. This convection is also extremely vigorous, and suspected to be one of the causes of mass-loss in RSGs. 
It must thus be understood in details. Evidence has been accumulated that there are asymmetries in the photospheres of RSGs, but 
detailed studies of granulation are still lacking. Interferometric observations offer an exciting possibility to tackle this question, but 
they are still often interpreted using smooth symmetrical limb-darkened intensity distributions, or very simple spotted ad hoc models. 
Aims. We explore the impact of the granulation on visibility curves and closure phases using the radiative transfer code OPTIM3D. 
We simultaneously assess how 3D simulations of convection in RSG with C0 5 BOLD can be tested against these observations. 
Methods. We use 3D radiative-hydrodynamics (RHD) simulations of convection to compute intensity maps at various wavelengths 
and time, from which we derive interferometric visibility amplitudes and phases. We study their behaviour with time, position angle, 
and wavelength, and compare them to observations of the RSG a Ori. 

Results. We provide average limb-darkening coefficients for RSGs. We detail the prospects for the detection and characterization 
of granulation (contrast, size) on RSGs. We demonstrate that our RHD simulations provide an excellent fit to existing interferomet- 
ric observation of a Ori, contrary to limb darkened disks. This confirms the existence of large convective cells on the surface of 
Betelgeuse. 

Key words, stars: supergiants - stars: atmospheres - hydrodynamics - radiative transfer - techniques: interferometric 

1 . Introduction l2007h . Although these MARCS models dGustafsson et alJl2008h 

give a good fit of the optical spectra allowing the derivation of 
Massive stars with masses between roughly 10 and 25 M Q spend T ^ and re ddening, problems remain. There is a blue-UV ex- 
some time as red supergiant (RSG) stars being the largest stars cess in many of the observed spectra, that may point to scatter- 
in the universe. They have effective temperatures, T eff , rang- ing by circumstellar dust, or to an insufficiency in the models, 
ing from 3 450 to 4 100 K, luminosities of 20 000 to 300 000 L xhere is also a m i smatc h in the IR colours, that could be due to 
and radii up to 1 500 R Q (|Levesque et al. || 2005|) . Their luminosi- atmosp heric temperature inho mogeneities characteristic of con- 
ties place them among the brightest stars, visible to very large vec tion dLevesque et al.ll2006h . 
distances. There is however a number of open issues. They 
shed large amounts of mass back to the interstellar medium, 
but their mass-loss mechanism is unidentifi ed, although Alfven 



In recent years, hydrodynamical modeling of convection in 



and acoustic waves have be en prop osed (lHartmann & Avrettl 
11984 !Pijpers&Hean3ll989t ICuntzl Il997h . as well as acous- 
tic wa ves and radiation pressure on molecules dJosselin & Plezl RSGs has lagged behind that of solar type stars du e to the neces- 
12007b . Their ch emical compositio n is l argely unknown, de spite sity to include the whole star in the simulation box.|Freytag eTaj] 
the work of e.g. ICarr et al.1 (l2000h . and lCunha etaD (120071) . due (|2Q02|) have succeeded in doing such numerical simulations of a 
to difficulties in analysing their spectra with broad, asymmet- typical RSG. We have thus engaged in an important effort to im- 
ric lines with variations suspected to stem from a convection prove our understanding and description of RSGs using detailed 
patt ern consisting of large gran ules and (super-)sonic veloci- numerical simulations and a large set of observational material, 
ties ( Josselin & Plez 2007; Grav 2008). Pro gress has been made 
recently, with their T^- scale being revised both at solar and 
Magellanic Clouds metallicities using ID hydrostatic models 

dLevesque et alJ |2^M|Masse^^|20a|Li7 e S que et al. | This paper is the first in tMs serfes and it is aimed tQ explore 

the granulation pattern of RSGs and its impact on interferometric 

Send offprint requests to: A. Chiavassa observations. 
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2. 3D radiative transfer in Radiative-hydrodynamics 
simulation 

2.1. 3D hydrodynamical simulations with C&BOLD 

The numerical simulatio ns employed in this work have been ob- 
tained using C0 5 BOLD dFrevtag et alJl2002h iFrevtag & Hofnel 
I2QQ8I) in the star-in-a-box configuration: the computational do- 
main is a cube, and the grid is equidistant in all directions. 
All six faces of the cube use the same open boundary con- 
ditions for material flows and emergent radiation. In addition, 
there is an "inner boundary condition": in a small spherical re- 
gion in the center of the cube a source term to the internal 
energy provides the stellar luminosity and a drag force brakes 
dipolar flows through it. Otherwise, the hydrodynamics and 
the radiative transfer scheme do not notice the core and inte- 
grate right through it. Radiation transport is strictly in LTE. 
The grey Rosseland mean opacity is a function of gas pres- 
sure and temperature. The necessary values are found by inter- 
polation in a table which has bee n merged at around 12000K 
from high-temperature OPAL data (llglesias et a 3l992h and low- 
temperature PHOENIX data (lHauschildt et alJll997b by Hans- 
G unter Ludwig. Some mor e technical information can be found 
in lFrevtag & Hofnerl(l2008l) . the C0 5 BOLD Online User Manual 
(www.astro.uu.se/~bf/co5bold_main.html), and in a forthcoming 
paper by Freytag (2009). 

The 12 M© model we use in this paper (st35gm03n07) is the 
result of intensive calculations which have led to about 7.5 years 
of simulated stellar time. It has a numerical resolution of 8.6 R© 
in a cube of 235 3 grid points. The model parameters are a lu- 
minosity of L = 93000 ± 1300 L , an effective temperature of 
T eff = 3490+ 13 K, a radius of R = 832+0.7 R , and followingly 
a surface gravity log g = -0.337 + 0.001. These values are aver- 
ages over spherical shells and over time (over the last year), and 
the errors are one sigma fluctuations with respect to the average 
over time. We define the stellar radius, R, and the effective tem- 
perature, r e ff, as follows. First, we compute the average tempera- 
ture and luminosity over spherical shells, T(r), and L(r). We then 
search the radius R for which L(R)/(47tR 2 ) = crT 4 (R), where cr 
is the Stefan-Boltzmann constant. The effective temperature is 
then T e ff = T(R). Fig. [T] shows the value of the radius, tempera- 
ture and luminosity over the last 3.5 years. The radius drifts by 
about -0.5% per year and seems to have stabilised to R = 832 R© 
in the last year. T Q s fluctuates by + 1 % over the whole sequence, 
with a constant average. The luminosity fluctuations are of the 
order of +4%, reflecting the temperature variations, with a de- 
crease of about 1% per year in the first years, reflecting the ra- 
dius decrease. These drifts indicate that the simulation has not 
completely converged in the first years. In this work we consider 
the whole 3.5 year sequence, despite the small radius drift, in 
order to have better statistics. The preceding 4 years of the simu- 
lation are not considered here, since they show larger drifts. The 
interferometric observables derived in this work are not sensitive 
to the drift of the parameters. 

This is our "best" RHD simulation so far because it has stel- 
lar parameters close st to real RSGs (e.g., 3650 K for a Ori, 
iLevesque et al.ll200"5l) . New simulations with different stellar pa- 
rameters are in progress and they will be analyzed in a forthcom- 
ing paper. 



2.2. Radiative transfer code: OPTIM3D 

We have developed a 3D pure LTE radiative transfer code, 
OPTIM3D to generate synthetic spectra and intensity maps 
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Fig. 1. Radius (Panel A), luminosity (Panel B) and temperature 
(Panel C) as a function of time for the simulation used in this 
work. The radius is fitted with the law: R(t) = 936.9 - 4.3 x t for 
t < 23.8 yrs years and R = 832 R© for t > 23.8 yrs. 



from snapshots of the 3D hydrodynamical simulations, taking 
into account the Doppler shifts caused by the convective mo- 
tions. The radiation transfer is calculated in detail using pre- 
tabul ated extinction coeffici ents generated with the MARCS 
code dGustafsson et al. 2008). These tables are functions of tem- 
perature, density and wavelength, and wer e computed with the 
solar composition of lAsplund et al.l (120061) . The tables include 
the same extensive atomic and molecular data as the MARCS 
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models. They were constructed with no micro-turbulence broad- 
ening and the temperature and density distribution is optimized 
to cover the values encountered in the outer layers of the RHD 
simulations. The wavelength resolution is R = A/AA = 500000 
and we checked that this resolution is sufficient to ensure an ac- 
curate calculation of broadened line profiles of RSGs even after 
interpolation of the opacity at the Doppler shifted wavelengths. 

The monochromatic intensity emerging towards the observer 
at a given position on the simulation can be computed by inte- 
grating the source function along a ray perpendicular to a face of 
the cube, at that position. In LTE it reads: 

h(0)= [ A S A {tA)e~ h dt A (1) 
Jo 

where Ix is the intensity, tx is the optical depth along the ray in- 
creasing inwards, Tx is the maximum optical depth reached along 
the line-of-sight, and S x = B A (T), the Planck function at the 
temperature T, is the source function. A Gauss-Laguerre quadra- 
ture of order n can be performed to evaluate the integral, Eq. (Q]), 
when Tx —> co. This method is much faster than a detailed inte- 
gration along the discretized ray, as it uses only the value of the 
source function at n depth points weighted with n predetermined 
weights. This method is reliable as long as the source function 
is sufficiently smooth along the optical depth scale, and is well 
known at the quadrature points. This is not always the case in 
our simulations where the optical depth scale may jump by large 
amounts between 2 successive cells, e.g. from r = 1 to r = 300 
for extreme cases. The source function must then be interpolated 
at intermediate optical depths, and the result is largely dependent 
on the way this interpolation is performed. Note however that 
this is also the case for a detailed integration, where too large 
jumps in the source function or optical depth scale will cause un- 
certainties in the resulting intensity. We checked for differences 
between a Gauss-Laguerre quadrature and a detailed summation 
of the contributions from all cells, with different kinds of inter- 
polations, and found differences in intensities emerging from a 
single ray of less than 10% on average, with some extreme cases 
reaching more than 100% due to a particularly illconditioning of 
the source function. The average differences being in an accept- 
able range, we therefore rely on the Gauss-Laguerre quadrature, 
with a linear interpolation of the source function on the logarith- 
mic r-sca le. The quadrature points and weights we use are listed 
in Tab. Q] ( Abramowitz & Stegu For the rays where the 

optical depth does not reach high enough values, we carry a de- 
tailed summation of the contribution from all cells along the ray. 

In practice, once the input simulation is read, OPTIM3D in- 
terpolates the opacity tables in temperature and logarithmic den- 
sity for all the simulation grid points using a bi-linear interpo- 
lation. The interpolation coefficients are computed only once, 
and stored. Bi-linear interpolation has been preferred to spline 
interpolation because: (i) spline is significantly more time con- 
suming, (ii) and comparisons with other codes do not show great 
improvements using splines (see below). Then, the logarithmic 
extinction coefficient is linearly interpolated at each Doppler- 
shifted wavelength in each cell along the ray, and the optical 
depth scale along the ray is calculated. Eq. (Q} is then integrated, 
giving the intensity emerging towards the observer at that wave- 
length and position. This calculation is performed for every line- 
of-sight perpendicular to the face of the computational box, and 
for all the required wavelengths. 

Comparisons with existing codes were carried out. The 
spectral synthesis c ode Turbospectrum (Plezetal ] fl993l 
I Alvarez & Plezj[T998[ and further improvements by Plez) was 
used with one-dimensional MARCS models, where the source 



Table 1. Gauss-Laguerre quadrature weights for n=10. 
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0.137793470540 


3.08441 115765^ 
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-04 


11.843785837900 


2.82592334960^ 


-05 


16.279257831378 


4.2493 1398496£ 


-07 


21.996585811981 


1.83956482398^ 


-09 


29.920697012274 


9.91182721961^ 


- 13 



fonction is very well sampled on the r-scale. OPTIM3D 
computations made with bi-linear interpolation deviate by 
less than 5%, and the deviation decreases to 0.2% with 
spline int erpolation. We als o checked OPTIM3D against 
Linfor3D ( Cayrel et al] 120071 for the Non-LTE version, and 
http://www.aip.de/~ mst/Linfor3D/linfor_3D_manual.pdf for the 
LTE version) using 3D C0 5 BOLD local models. We compared 
synthetic spectra computed for three artificial iron lines (with 
increasing strength) centered at a laboratory wavelength of 5500 
A and using the same abundances. The discrepancy between the 
results of the codes is less than 3% and it becomes even less 
than 0.2% when a spline interpolation of the opacity tables is 
used in OPTIM3D (with a significant increase of the CPU time). 
Finally, comparis ons were made w ith the spectral line formation 
code used by, e. g. jAsplundl (l200Qh for 3D local convection si mu- 
lations carried out with the code bv lStein & Nordlundl(ll998h for 
giant stars (ICollet et alJ2007h . The tests have been carried out on 
the [OI] line at 6300.3 A and various Fe I and Fe II lines around 
5000 A. The discrepancies between the resulting synthetic spec- 
tra are less than 2%, and become even less than 0.6% when a 
spline interpolation of the opacity tables is used in OPTIM3D. 
Thus, the interpolation is the main source of error. In conclusion, 
if only a few lines are computed for, e.g., accurate abundance de- 
terminations, Linfor3D or the Asplund code are better because 
they mostly avoid interpolations into opacity tables. The code 
used by Asplund performs bi-cubic interpolations of the contin- 
uum opacity and of the individual number densities, whereas we 
interpolate the total opacity from all lines contributing at a given 
wavelength, which is in principle less accurate. So, when a large 
wavelength range must be calculated taking into account many 
molecular and atomic lines simultaneously, OPTIM3D is a bet- 
ter, faster choice, that still provides a result accurate to a few 
percents. 

3. Simulated images in the H and K bands: giant 
convective cells 

In this work, we analyze the properties of the simulations in the 
H and K bands where many interferometric observations have 
been done, and existing interferometers, e.g. VLTI/AMBER, 
operate routinely. We calculated intensity maps for a series of 
snapshots about 23 days apart covering 3.5 years of the model 
described above. We use the transmission curve of the four K 
band filters mou nted on FLUOR (Fiber Linked Unit for Optical 
Recombination; Coude Du Fore sto et al. and of the H 

band filter mounted on I ONIC (Integrated O ptics Near-infrared 
Interferometric Camera; Berger et all 120031) at the IOTA inter- 
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Fig. 2. The transmission curves of the 4 narrow band filters 
mounted on the FLUOR instrument at IOTA together with the K 
band synthetic spectrum of a snapshot of the simulation and the 
corresponding continuum (bottom black curve). From the top, 
the spectra computed with only H 2 (red), only CO (green), 
and only CN (blue) are shown with an offset of respectively 
0.9 x 10 33 , 0.6 x 10 33 and 0.3 x 10 33 erg cm" 2 s" 1 A" 1 . 
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Fig. 3. The transmission curve of the filter mounted on IONIC at 
IOTA together with the H band synthetic spectra computed as in 
Fig. [2l From the top, the offset of the spectra is 3 x 10 33 , 2 x 10 33 , 
and 1 x 10 33 erg cm" 2 s" 1 A" 1 . 



ferometer dTraub et al.ll2003h . The K band filters (Fig. 0) are: 
K203 (with a central wavelength of 2.03 /mi), K215 (2.15 /mi), 
K222 (2.2 /mi), and K239 (2.39 /mi). The H band filter has a 
central wavelength of 1.64 /mi (Fig. [3]). The resulting intensities 
reported in this work are normalized to the filter transmission as: 

J I A T(A)dA 



f T(A)dA 



where h is the intensity and T (A) is the transmission 



curve of the filter at a certain wavelength. The intensity maps are 
showed after applying a median [3x3] smoothing (see Section©. 

It can be seen from our simulations (see Fig. [4]) that the 
surface of the stellar model is covered by few large convective 
cells of a size of about 400 to 500 R© that evolve on a time- 
scale of years. These cells have strong downdrafts that can pen- 
etrate down to the stellar core (Freytag et al. 2002, and 2009, 
in preparation). Near the surface, there are short-lived (a few 



months to one year) s mall-scale (50 to 10 R©) granules (bot- 
tom panels of Fig. 0]). iFreytag et al] (1 19971) found a relation be- 
tween the mean horizontal size of convective granules x gran and 
the atmospheric pressure scale-height defined as H p o = for 
GK dwarfs and subgiants. It is unclear if such a relation can 
be extrapolated to 3D simulations of RSGs. Using it we find 
%an/^* = 10 x H p o/R+ = 0.1, for parameters appropriate for 
a RSG atmosphere dominated by gas pressure. Obviously, this 
leads to a size muc h smaller than what can be seen in Fig. [4] 
IFreytag et al.l (1 19971) found that a value of x gmn /H p o = 10 would 
fit 2D simulations for GK dwarfs and subgiants, but they show 
also that A-type and F-type stars lie above the curve indicating 
that they have larger granules. These stars have large turbulent 
pressure that may dominate over the gas pressure in turbulent 
convective layers. Following Gustafss on et al.l (120081) . we write 
Aurb = fipVt mh , where v tur b is the turbulent velocity, p is the gas 
density, and j3 is a parameter close to one, whose value depends 
on the anisotropy of the velocity field. A better way to express 

H p0 is thus Hpo = g£ (l +/3y where y is the adia- 

batic exponent, and c s the sound speed. If v tur b is only a factor 2 
larger than c s , H p o is increased by a factor of about 5. This is the 
case for our RSG simulation where /\urb/^gas ~ 2 at the surface, 
R*, as deter mined in Sect . 12. II T his gives then x gra n/^*=0.5, ex- 
trapolating IFreytag et al] (1 19971) formula. This is more consis- 
tent with the large granules visible on intensity maps in Fig. 01 
There are further mechanisms that might influence the size of the 
granules: (i) in RSGs, most of the downdrafts will not grow fast 
enough to reach any significant depth before they are swept into 
the existing deep and strong downdrafts enhancing the strength 
of neighboring downdrafts; (ii) radiative effects and smoothing 
of small fluctuations can cause an enhancement of growth time 
for small downdrafts while the granule crossing time is short due 
to lar ge horizontal velocities; (ii i) sp hericity effects, see for ex - 
am plelFrevtag & Ludwid (l2007l) . and lSteffen & Frevtagl (120071) : 
(iv) IFreytag et al] (119971) use the effective temperature and the 
pressure scale height at the bottom of the photosphere as refer- 
ence, however, also layers below the photosphere can matter; (v) 
numerical resolution (or lack of it) could matter. 



4. Intensity profiles 

The simulated RSG atmospheres appear very irregular, perme- 
ated with structure and dynamics. The surface inhomogeneities 
and their temporal evolution induce strong variations of the 
emerging spectra, and intensity profiles. In this section we anal- 
yse the average centre-to-limb intensity profiles, and their time 
variations. 



4.1. Surface inhomogeneities and temporal evolution 

The top left panel of Fig. [5] shows a three-dimensional image rep- 
resentation of the intensity emerging from one face of a snapshot 
of the simulation in the H band. The K band appearance is sim- 
ilar. This image shows very sharp intensity peaks two to three 
pixels wide. This is also noticeable in the top right panel of the 
Figure as small bright (up to 40% brighter than the surrounding 
points) patches. These patches result from the ill-conditioning of 
the source function, due to the lack of spatial resolution around 
Tx = 1 along some lines of sight where the source function may 
have a large jump (see Sect. 12.2b . Attempts have been made to 
solve this problem through the interpolation of the source func- 
tion and opacity inside CO 5 BOLD, but they caused numerical 
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Fig. 4. Top 6 panels: maps of the intensity in the IONIC filter (linear scale with a range of [0;2.5 x 10 5 ] erg cm -2 s" 1 A -1 ). The 
different panels correspond to snapshots separated by 230 days (~ 3.5 years covered). Bottom 6 panels: successive snapshots 
separated by 23 days (~ 140 days covered). 
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instabilities. The unique solution is to increase the number of 
grid points, and that necessitates larger and faster computers. 

Radial intensity profiles within a given snapshot show large 
variations with position angle of their radial extension of about 
10% (see bottom left panel of Fig.0. The variation with time of 
the intensity profiles are of the same order of magnitude (10%, 
see bottom right panel of Figure). 

4.2. The limb darkening law 

Despite the large azimuthal variations of the intensity profiles, 
and their temporal variations, it is interesting to derive radially 
averaged intensity profiles for each snapshot. These may be be 
used, e.g., as a first approximation to interpret interferometric 
observations, in replacement o f limb-darken ing (LD) laws com- 
puted from hydrostatic models (IClaretll2000l) . Bottom right panel 
of Fig. [5] shows all the radially averaged intensity profiles ob- 
tained from the simulation. 

We use a LD law of the form : 

where is the intensity, are the LD coefficients and \i = 
cosO with the angle between the line of sight and the radial 
direction, ji is related to the impact parameter r/R* through 
r/R+ = y/l - jj 2 , where R* is the stellar radius determined as 
in Sect. 12.11 The average intensity profiles were constructed us- 
ing rings regularly spaced in fi for \i < 1 (i.e., r/R+ < 1), and 
adding a few points for r > R+ up to the numerical box limit. The 
standard deviation of the average intensity, 07^), was computed 
within each ring. There is a small tail at r > R+ that gives a mi- 
nor contribution to the total flux (less than 1%, see bottom right 
panel of Fig. [5]), and cannot be fitted with Eq.O We fitted the ra- 
dially average profiles of all the snapshots of the simulation (57 
profiles 23 days apart covering 3.5 years). The fit was weighted 
by 1/07^) in order to decrease the importance of central points 
with poor statistics. The fit was first made on IQi)/I noTm , with 
^norm the intensity normalized to the area subtended by the curve: 
/norm = 7 ( r /^*) j^[ s was ^one m or d er t diminish the 

£ Kr/RJdr/R* 

impact of intensity fluctuations between snapshots on the fitting 
coefficients. 

In Tab. |2 we give the values of the four LD coefficients aver- 
aged over all 3.5 years, and renormalized to disk center, for the 
IONIC H-band filter, and for the K222 filter (because the sensi- 
tivity of the FLUOR instru ment is always bette r in the contin- 
uum than in molecular bands IPerrin et al.ll2004bl and it samples 
the maximum transmission region of the K band). Fig. [6] shows 
an example of LD fit. The intensity profiles for different position 
angles for the same snapshot being very different (see Fig. [5j 
bottom left panel), the fitting coefficients are very scattered. The 
time averaged LD fits give however an indication of the shape of 
the intensity profile in the H and K bands (note that they are very 
similar). They are of course very different from simple first or- 
der LP laws. They also differ from LD laws calculated by lClaretl 
(120001) for parameters appropriate for RSGs (see Fig. [7]). When 
comparing to observations of RSGs, we recommend to use our 
fits. Ideally one should use single snapshots as we do below in 
our analysis of Betelgeuse, as they may deviate from the average 
LD fit by large amounts (see Tab. [2]). 
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Fig. 6. Example of a LD fit (dashed line) using the LD law 
described in the text for the radially averaged intensity profile 
(solid line) emphasized in Fig. [5] (bottom right panel). The in- 
tensity is normalized to the area subtended by the curve. This 
best fit has a x 2 = 0.02. 
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Fig. 7. The time averaged H-band radial intensity profile of our 
simulation (solid line), and the fit of Tab. [2] (dashed line). A fully 
LD (dash-dotted li ne), a partiall y LD (triple dot-dashed line), 
and a LD fit from IClareT (I200QL dotted line) for comparable 
RSG parameters are plotted for comparison. 



Table 2. Time-averaged limb-darkening coefficients for the 
RHD simulation (see Eq.0 cr is the standard deviation over time 
(57 profiles covering 3.5 years. 
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a central wavelength of the corresponding IONIC filter 

b central wavelength of the corresponding K222-FLUOR filter 
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Fig. 5. Top left panel: three-dimensional image of a snapshot from Fig. HI Top right panel: Intensity map of the same snapshot 
represented using the histogram equalization algorithm in order to underline the thin bright patches due to the undersampling of the 
r scale. Bottom left panel: intensity profiles for three position angles of the same snapshot. The numerical box edge is at impact 
parameter rlR+ ~ 1.3. The intensity is normalized to the intensity at disk center. Bottom right panel: radially averaged intensity 
profiles for all the snapshots of the simulation (grey); one snapshot of the simulation is emphasized with a solid black line. The 

rl 3 

intensity is normalized to the area subtended by the curve, area= J Q / (r/R+) drjR*. 



5. Visibility curves and phases 

5.1. Computation 

The granulation pattern has a significant impact on interferomet- 
ric visibility curves and phases. We try here to derive their char- 
acteristic signature. 

We compute visibilities and phases using the IDL data vi- 
sualization and analysis platform. For each image, a discrete 
Fourier transform is calculated. In order to reduce the prob- 
lems due to the finite size of the object and to avoid edge ef- 
fects, the resolution in the Fourier plane is increased by padding 
the input 235 x 235 pixels image with zeros up to a size of 
2048 x 2048 pixels. The visibility V is defined as the modu- 
lus \z\, of the complex Fourier transform, z - x + iy, normalized 
to the modulus at the origin of the frequency plane, \zo\, with 
the phase tan# = 3(z)/2l(z). When dealing with observations, 



the natural spatial frequency unit is arcsec . As we study theo- 
retical models, we use instead R" 1 units. The conversion factor 
between these is 

V [arcsec~ l ] = V [R" 1 ] • d [pc] • 214.9, (3) 

where 214.9 is the astronomical unit expressed in solar radii, 
and d is the distance of the observed star. The relation between 
the baseline, B (in m), of an interferometer, and the spatial fre- 
quency v (in arcsec -1 ) observed at a wavelength A (in /dm) is 
v = 5//1/0.206265. As our calculated images are affected by the 
source function jumps, we investigated how the visibility curves 
are affected by the resulting bright spikes (Fig. [5]). We compare 
in Fig. [8] the visibility curves computed for one projected base- 
line from the raw image, and after applying a median [3x3] 
smoothing effectively erasing the artifacts. The visibility curves 
are affected by these artifacts mostly for frequencies greater than 
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0.03 R" 1 (corresponding to 33 R©, i.e., ~4 pixels). We can there- 
fore apply this cosmetic median filter, as it will not affect the 
visibilities at lower frequencies, that are the only ones to be ob- 
servable in practice with modern interferometers. 

1000 



1.0000 




0.0001 1 ■ i\...A\i.i.A^..[i.\im.m..i 

0.00 0.01 0.02 0.03 0.04 0.05 
Re" 1 



Fig. 8. The solid line is the visibility curve for the IONIC filter 
intensity map of Fig. [5] (top right). The dotted line is computed 
for the same map after applying a [3 x 3] median smoothing. 

We now study the first few lobes of the visibility curves of 
our simulations, and how they are affected by asymmetries and 
surface structure. 

5.2. The first lobe 

The first lobe of the visibility curve is mostly sensitive to the 
radial extension of the observed source. Fig. [9] (bottom panel) 
shows the visibility curves computed for 36 different angles from 
the intensity map of one snapshot in the IONIC filter (top panel). 
A dispersion of the visibility curves (thin grey lines in Fig. [9]) is 
noticeable. This behavior is similar for all the snapshots. These 
synthetic visibilities have been compared to a uniform disk (UD) 
model (solid line in Fig. O, and with limb-darkened (LD) mod- 
els. We use both a fully limb-darkened disk (I^/h = yu, dotted- 
dashed line in the Figure), and a partially limb-darkened model 
with a\ - -0.5 (I/j/h = 0.5 + 0.5 • \i, dashed line in the Figure). 
The radius determined by fitting a UD disk model to the com- 
puted visibilities ranges from 794 to 845 R© for the 36 angles, 
up to 5% smaller than R+ =836.5 R©, the radius of the simulation 
determined as described in Sect. I2.U The partially-, and fully- 
darkened models are respectively ~ 2%, and only -1% smaller 
than R+. In Fig. [9] is also shown the visibility amplitude result- 
ing from our average LD fit of Tab. [2] The resulting diameter 
is then 842 Rq very close to the simulation radius. Note, that 
iNardetto et~aiT (120061) also found that the UD radius is about 4 
to 5% smaller than the photospheric radius of their simulation 
of Cepheids, and that the LD radius is much closer to the radius 
of their simulation. Stellar diameters determined with UD or LD 
fits of observed first visibility lobe of RSGs will be affected by 
these systematic errors. As will be shown below, observations of 
higher spatial frequencies will greatly improve the knowledge 
of the limb-darkening, and of asymetries, thus helping in better 
constraining the radius as well. 

It is interesting to compare the angular and temporal visibil- 
ity fluctuations at one sigma, defined as F = cr/ vis: (i) the tem- 
poral evolution, fixing one angle and following the RHD simu- 
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Fig. 9. Top panel: intensity map in the IONIC filter (the range is 
[0;2.5 x 10 5 ] erg cm -2 s" 1 A -1 ). Bottom panel: visibility curves 
from the above snapshot computed for 36 different angles 5° 
apart (thin grey lines). Note the logarithm visibility scale. The 
solid black curve is a UD model, with a radius of 810 R . 
The dashed black line is a partially LD disk with a radius of 
822 R©. The dot-dashed line is a fully LD disk with a radius 
of 830 R . the triple-dot-dashed line is our average LD law (cf. 
Tab. [2]) for a radius of 842 R©. The stellar parameters of this 
snapshot are: L = 98 4OOL , R* = 836.5 R Q , r eff = 3534 K and 
log g = -0.34. 



lation for 3.5 years with a time-step of ~23 days; (ii) the angular 
evolution, considering a single snapshot and computing the visi- 
bilities for 36 different angles 5° apart. Fig. [TO] shows that, in the 
first lobe, temporal and angular fluctuations have the same or- 
der of magnitude. The fluctuations are less than 1% at frequency 
-0.00040 Rq 1 (at this frequency, the visibility is greater than 
50%), they are ~ 3% at frequency 0.00057 R" 1 , and are close to 
~ 10% at 0.00069 R" 1 . 
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Fig. 10. standard deviation of the visibility normalized to the vis- 
ibility in the first lobe. The solid line indicates the temporal fluc- 
tuations for one fixed angle over 3.5 years. The trend is similar 
for the other angles. The dashed curve corresponds to the angular 
fluctuations of the snapshot in Fig. [9] 



0.14 - 




0.0010 0.0015 0.0020 0.0025 
Re" 1 

Fig. 11. Same as Fig. [9] for the second, third, and fourth lobes. 
In addition, the dotted line is the visibility curve for the particu- 
lar azimuth parallel to the x-axis of the IONIC intensity map of 
Fig.H 



5.3. The second, third and fourth lobes: signature of the 
convection 

As in Sect. 15.21 we analyze the angular, and temporal visibility 
fluctuations at one sigma with respect to the average value in the 
H band (IONIC filter). Fig. [TT] shows an enlargement of the the 
second, third and fourth lobes of the visibility curves computed 
for different position angles. The dispersion increases clearly 
with spatial frequency, and visibilities deviate greatly from the 
UD or LD cases, due to the small scale structure on the model 
stellar disk. The same is true for temporal fluctuations of the vis- 
ibility at a given position angle. Fig. [12] shows the temporal fluc- 
tuations of the visibilities for one fixed position angle, as well as 
the angular fluctuations for the snapshot of Fig.[9j As for the first 
lobe, there is no clear distinction between angular and tempo- 
ral fluctuations. Relative fluctuations are of course large around 
the minima of visibility, where visibilities are also more difficult 
to measure. However, with the precision of current interferome- 
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Fig. 12. Same as Fig.[l0lfor the second, third, and fourth lobe. 



ters (e.g., 1% for visibilities of ~ 5 - 10% for VLTI- AMBER), 
it should be possible to characterize the granulation pattern on 
RSGs. This requires observing the third and the fourth lobes and 
not limiting the observation at the first and second lobes, that 
give only an information on the radius and LD. The signal to be 
expected in these lobes is higher than the UD or LD predictions 
(see dashed line in Fig. [TT]). Efforts should therefore be put on 
observing at these frequencies. 

It may however turn out that approximations in our mod- 
elling (e.g., limited spatial resolution, grey radiative trans- 
fer) significantly affect the intensity contrast of the granula- 
tion. Indeed, the radiation transfer in our RHD models uses a 
frequency-independent grey treatment to speed up the calcula- 
tions. This approximation leads to errors in the mean temper- 
ature structure in the optically thin layers that are difficult to 
quantify. The implementation of non-grey opacities can decrease 
the temperature flu ctuations compared to the grey case (e.g., 
iLudwig et al.ll 1994 for local RHD models). As a consequence, 
the intensity contrast will be decreased, reducing the visibility 
fluctuations. To investigate its impact on visibilities, we artifi- 
cially decrease the intensity contrast on one of our images, we 
use the snapshot with its nominal intensities as reference. We 
first fit a LD law (as in Sect [4]) to the radially average inten- 
sity profile. After subtracting this average profile from the in- 
tensity map we are left with the fluctuations caused by granu- 
lation. We measure the contrast C re f= ! max ~! min . It is then easy to 

-•max """■'min 

scale that contrast before adding again the LD profile, to recover 
a reduced contrast image. An example of the resulting intensi- 
ties is shown in Fig. [13] (top row). At a contrast of only 1% of 
the nominal one, small surface structures are hardly visible. As 
previously, we determine <x v i S /vis for all the images with vari- 
ous contrats, around the top of the second (-0.0010 R© 1 ), third 
(-0.0016 R" 1 ), and fourth lobes (-0.0022 R" 1 ). The bottom left 
panel of Fig. [13] shows that when the contrast is reduced, and 
the surface structures fade out, the resulting visibility fluctua- 
tions decrease similarly in all the lobes (almost proportionally 
to the intensity contrast decrease). Reducing the contrast brings 
of course the visibility curves towards the visibility of the LD 
profile (Fig. [13] bottom right panel). This proportionality can be 
used to determine the granulation contrast from observations of 
the visibility fluctuations with time or position angle. 
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Fig. 13. Top left panel: three-dimensional image of a snapshot with nominal intensities. Top right panel: same snapshot with a 
feature contrast reduced to 1%. Bottom left panel: standard deviation of the visibility curves at 36 angles 5°) apart, for decreasing 
feature contrast. Solid line is for the top of the second lobe (~ 0.0010 R© 1 ), dashed line is for the top of the third lobe (~ 0.0016 R© 1 ), 
and dotted line is for the top of the fourth lobe (~ 0.0022 R© 1 ). Bottom right panel: Visibility in the second, third, and fourth lobes 
for one particular position angle. The dot-dashed line shows the original simulation contrast. The dashed line, and the dotted line 
show the visibility with a feature contrast reduced to 50%, and 1% respectively. The solid line is the fitted LD profile computed for 
this snapshot. 



5.4. Importance of spectral resolution in interferometry: the 
H and K bands 

Interferometric observations done through a broad band filter 
blend information from the lines and continuum. Spectral res- 
olution allows to recover much richer information, both from 
visibility moduli and phases. The VLTI- AMBER interferome- 
ter provides spectral resolutions of R=35, 1500, and 12000. 
In order to show the differences between these resolutions, we 
compute intensity maps around the CO first overtone line at 
23041.75 A (log(g/)=-5.527 and ^=0.180 eV) for the three 
resolutions (Fig.[T4l). The resulting images are shown in the cen- 
tral row of the Figure, and the spectrum in the top row. The con- 



trast, defined as in previous Section, is similar for the low and 
medium resolution images but it is ~ 30% lower in the high res- 
olution image, at the CO line wavelength. Large fluctuations are 
seen in all lobes, but they are smaller than those seen in the H 
band IONIC filter (see Fig [12]). The visibility fluctuations for the 
high spectral resolution image in the CO line are larger (second, 
third and fourth lobes; see dotted line in bottom panel), despite 
a lower intensity contrast, presumably because of the darkening 
of large patches of the simulated stellar surface. 

We also computed wavelength dependent visibility curves in 
the H band for the high and medium VLTI- AMBER resolutions. 
Fig. [I5]displays a three-dimensional view of the visibility curves 
with a resolution of 12000, and 1 500 (top panels). The simu- 
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Fig. 14. Top row: synthetic spectrum centered on the CO first overtone band head. The left panel shows the range of wavelengths 
spanned by one resolution element at the VLTI- AMBER low spectral resolution of 35. The right panel shows the same for the 
VLTI- AMBER medium spectral resolution of 1 500, and the high spectral resolution of 12 000 (thick mark). Central row: intensity 
maps for those three spectral resolution elements. The intensity range is [0;10 5 ] erg cm -2 s" 1 A -1 . Bottom row: standard deviation of 
the visibility in the second, third and fourth lobes Solid and dashed lines correspond low and medium resolution respectively, and 
the dotted line to high resolution. 
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lated star has been scaled to an a pparent diameter of -43.6 mas 
(the observed diameter of a Ori P errin et al. The dis- 

placement of the zero points with wavelength is easily seen, as 
well as the amplitude variations in the higher frequency lobes. 

In order to mimic differential observations with an interfer- 
ometer at medium and high spectral resolution, we also show 
in Fig. Q21 the variation of the visibility modulus with wave- 
length for a fixed baseline (15m, i.e. in the second lobe, at v = 
45arcsec _1 ). The visibility shows variations correlated with the 
flux spectrum: it decreases in absorption lines. In fact, at these 
wavelengths wiggles and dark spots appear on intensity maps 
(Fig. [T4] central right panel) increasing the visibility signal at 
frequencies higher than the second lobe. The visibility variations 
are much attenuated at lower spectral resolution. Observations 
at wavelengths in a spectral line, and in the nearby continuum 
will probe different atmospheric depths, and thus layers at differ- 
ent temperatures. They will thus provide important information 
on the wavelength dependence of limb darkening. Moreover, as 
the horizontal temperature and density fluctuations depend on 
the depth in the atmosphere, differential observations, with rela- 
tive phase determination will provide unique constraints on the 
granulation pattern. The visibility variations in Fig [15] such as 
the steep visibility jump from 0.123 to 0.107 between 1.5975 
and 1.5980 /mi, could be measured in differential interferomet- 
ric mode at high spectral resolution with the current precision 
at VLTI-AMBER (1% for visibilities of -5-10%), with optimal 
sky conditions. 

5.5. Closure phase: departure from circular symmetry 

As terrestrial atmospheric turbulence affects the phases of the 
complex visibilities with random errors, it is impossible to derive 
them for individual pairs of telescopes. Instead one uses closure 
phase between three telescopes, as the sum of all phase differ- 
ences removes the atmospheric contribution, leaving the phase 
inform ation of the object visibility untouched (see e.g jMonnierl 
|2007|) . The closure phase is thus an important complementary 
piece of information, which can reveal asymmetries of RSG at- 
mospheres. Fig. [16] shows the scatter plot of the closure phase 
of one snapshot of the RHD model computed in the IONIC fil- 
ter (the scatter is similar for the K222 filter). The behavior is 
similar for all the snapshots. We used 500 random baseline tri- 
angles with a maximum linear extension of 40 m, and plot the 
closure phase as a function of the triangle maximum baseline. 
The closure phases deviate from zero or ±n already at ~10 m 
(0.0008 R" 1 if we scale the model to an apparent diameter of 
43.6 mas at a distance of 174.3 pc). At higher baselines it is 
clearly different from zero or ±n, values which indicate a point 
symmetric brightness distribution. This is a clear signature of 
surface inhomogenities. The characteristic size distribution on 
the stellar surface can also be derived from the closure phase: the 
contribution of small scale convection-related surface structures 
increases with frequency. The first deviation at -0.0008 Rq 1 (just 
beyond the first zero, see Fig. [TT} corresponds to the deviation 
from circular symmetry of the stellar disk. It may be very effi- 
cient to constrain the level of asymmetry of RSG atmospheres by 
accumulating statistics on closure phase at short and long base- 
lines, as they are easily measured with great precision. Small de- 
parture from zero will immediately reveal departure from sym- 
metry. 

We also computed the closure phase for the different K band 
VLTI-AMBER spectral resolution intensity maps of Fig.[l4j The 
large deviations from circular symmetry are already noticeable 
at low spectral resolution (Fig. [T7] left panels) and the closure 
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Fig. 16. Scatter plot of closure phases in the IONIC filter cen- 
tered at 1.64 /im of 500 random baseline triangles with a maxi- 
mum linear extension of 40 m. Closure phases are plotted against 
the longest baseline of the triangle. The upper x-axis corresponds 
to synthetic observations of the simulation at an apparent di- 
ameter of 43.6 mas (which corresponds to a Ori at a distance 
of 174.3 pc). The axisymmetric case is represented by the grey 
lines. 



phase scatter do not differ much from the high spectral resolution 
one (right panel). This offers prospects of detecting asymmetries 
due to granulation without resorting to high spectral resolution. 



6. Comparison with interferometric observations of 
a Ori 

A first confrontation of our model predictions to real observa- 
tions is possible for a Ori. We compare the synthetic visibilities 
derived from our RHD simulations in t he continuum filter K 222 
(Fig. [2]) with the observation of a Ori by lPerrin et al.l (l2004al) that 
reach the third lobe in the K band. The absolute model dimen- 
sions have been scaled to match the interferometric observation 
in the first lobe. This corresponds to an apparent diameter of 
43.6 mas at a dis tance of 179 pc. These values are in agreement 
with Perrin e t all w ho found a diameter of 43.64 ±0.10 mas, and 
lHarper et~ai1 ([2008), who reported a distance of 197 ± 45 pc. 

We computed over 2000 visibility curves, and we find 
that the data are within the visibility fluctuations due to the 
granulation of t he sim ulation (Fig. [18]), as already shown in 
IChiavassa et al.l 120071 Within this large number of visibility 
curves, we find some that match all the observation points better 
than th e uniform disk (with a diameter of 43.33 mas: lPerrin et al.l 
l2004a|) . or limb-darkened disk model (linear limb darkening law, 
/ (jj) = 1 - a (1 - u ), with a diameter of 43.64 mas and a = 0.09 
also in lPerrin et al.h . See Fig. [181 The best match has a reduced 
X 2 =0.2l, and all the visibility curves fall within a x 2 range of 
[0.21,18.1]. Our RHD simulations are a great improvement over 
parametric models (the UD model with reduced x 2= 19.9, and 
the LD model with x 2 =22. 3) for the interpretation of these in- 
terferometric observations. The observations points in the first, 
second, and third lobes can be reproduced with a single visibility 
curve, from the projection at a particular position angle of one 
of our snapshots (see Fig. [18]). There is one observed point in 
the first lobe at 24.5 arcsec -1 which is difficult to reproduce. In 
fact, adjustments on the absolute model dimensions of the star in 
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Fig. 15. Top left panel: three-dimensional view of the visibility curves as a function of wavelength for a particular position angle. 
The spectral resolution is 12 000. Top right panel: same as in top left panel at a spectral resolution of 1 500. Bottom panel: Visibility 
as a function of wavelength for a baseline of ~15 m, i.e. the top of the second lobe, at one particular position angle. The simulation 
has been scaled to an apparent diameter of -43.6 mas. The synthetic spectrum convolved to a resolution of 12000 is over-plotted 
(thin solid line). The small black dots correspond to the highest resolution with AMBER (12 000), while the big red dots correspond 
to the medium resolution (1 500). The crosses show the uniform disk of 43.6 mas. When changing the position angle, the expected 
standard deviation is about 10% of the visibility (see Fig.fT2t. 



order to fit this point, would lead to mismatch of the other obser- 
vations at higher frequencies. However, this may be a problem 
with the calibration of the observation. 

A more detailed comparison with a Ori data in the H band 
(lHaubois et al.ll2006l) will be presented in a forthcoming paper 
(Chiavassa et al. 2009, in prep.) 

7. Conclusions 

Our radiation hydrodynamics simulations confirm that only a 
few large granules cover the surface of RSG stars. The gran- 
ules of the simulation we analyze here are 400-500 R in diam- 
eter, and have lifetimes of years. Smaller scale structures develop 
and evolve within these large granules, on shorter timescales (a 
month). 

We demonstrate that RHD simulations are necessary for a 
proper quantitative analysis of interferometric observations of 



the surface of RSGs beyond the smooth, symmetrical, limb- 
darkened intensity profiles. We give new average limb darkening 
coefficients within the H and K bands, that are significantly dif- 
ferent from commonly used UD or LD profiles. However, these 
LD coefficients fluctuate with time, and the average is only in- 
dicative. Our model surface granulation causes angular and tem- 
poral variations of visibility amplitudes and phases. In the first 
lobe, sensitive to the radius, fluctuations can be as high as 5%, 
and radii determinations can be affected to that extent: the radius 
determined with a UD fit is 3-5% smaller than the radius of the 
simulation, while the radius determined with a fully LD model 
is 1% smaller. 

The second, third, and fourth lobes, that carry the signature 
of limb-darkening, and of smaller scale structure, are very dif- 
ferent from the simple LD case. The visibility amplitudes can be 
greater than the UD or LD case, and closure phases largely differ 
from and ±n, due to the departure from circular symmetry. The 
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Fig. 17. Scatter plot of closure phases (cf Fig. [TS] for details) obtained from the VLTI- AMBER K band low, and high spectral 
resolution intensity maps (Fig.[l4l central left, and right panels). 



visibilities also show fluctuations with position angle, and with 
time, that are directly related to the granulation contrast. We also 
want to stress that high spectral resolution provides an extremely 
valuable information. The stellar surface appears dramatically 
different in an absorption line and in the nearby continuum, and 
differential observations should be easier to carry out with high 
precision. At lower resolution (e.g. R = 1 500) continuum and 
line information get mixed and there is a considerable loss of 
differential signal. 

We conclude that the detection of the signature of granula- 
tion, as predicted by our simulation, is measurable with todays 
interferometers, if observations of both amplitude and closure 
phase are made at high spatial frequencies (second, third, and 
fourth lobes, or even further if possible). These observations will 
provide direct information on the time scale of variation, and of 
the size and contrast of granulation. 

A few RSGs are prime targets for interferometry, thanks to 
their large diameter, proximity, and to a high infrared luminosity. 
Only 4 or 5 are sufficiently close and bright that imaging can be 
attempted, but a larger number (10-20) are within reach of less 
ambitious programs, like closure phase measurements. 

Three approaches can be combined to characterize the gran- 
ulation pattern: 

- searching for angular visibility variations, observing with the 
same telescope configuration (covering high spatial frequen- 
cies) and using the Earth rotation in order to span 6-7 posi- 
tion angles in one night should be enough, if measurement 
errors can be kept below 10%, for visibilities of the order of 
5 to 10%. One or two other telescope configurations would 
provide more frequency points, but then the change of con- 
figuration must be made within days, which is actually pos- 
sible at VLTI; 

- looking for temporal visibility fluctuations by observing at 
two (preferably more) epochs ~ 1 month apart with the same 
telescope configuration. This can be easily scheduled on ex- 
isting interferometer like CHARA, or the VLTI; 

- looking for visibilities as a function of wavelength, at high 
spectral resolution, in different spectral regions belonging 
to spectral features and continuum. If measurement errors 
can be kept close to 1%, for visibility of the order of -10%, 
variations of the visibility correlated with the flux spectrum 



could be detected, indicating variations of the radius, of the 
limb-darkening, or of the granulation pattern. Such relative 
measurements are more easily done at the required precision 
than absolute visibility measurements. 

These observations will bring us a wealth of information on 
the stars, but also on our RHD models. We know they suffer from 
limitations. The confrontation to observation will help us decide 
what approximations must be relaxed. The simulations are pri- 
marily constrained by execution time, which neccesitates sev- 
eral approximations, the most important ones being the limited 
spatial resolution and the complete lack of wavelength resolu- 
tion, i.e., grey radiative transfer. This speeds up the simulations 
to managable execution times of several months to one year of 
intensive calculations for about seven years of stellar time. A 
higher numerical resolution shows smaller scale structures ap- 
pearing within the g ranules that are already present in lower res- 
olution simulations (iChiavassa et al.ll2006l Fig. 10). This, how- 
ever, should not affect the first few visibility lobes. 

The approximation of grey radiative transfer is justified only 
in the stellar interior and it is not appropriate in the optically 
thin layers. The implementation of non-grey opacities (e.g., five 
wavelength groups employed to describe the wavelength depen- 
dency of the radiation field w ithin a multigroup radiative transfer 
scheme, see iNordlundll 19821 for details) would be an important 
improv ement for the hydrodynamical simulations. iLudwig et al.l 
(1 19941) found for local RHD simulations that frequency depen- 
dent radiative transfer causes an intensified heat exchange of a 
fluid element with its environment tending to reduce the temper- 
ature differences. Consequently, the temperature fluctuations in 
the non-grey local models are smaller than in the grey case. This 
is also expected for global RSG models, and the result of such 
a decrease of the temperature fluctuations, would be a decreased 
intensity contrast and decreased visibility fluctuations. 
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Fig. 18. Comparison of the RHD simulation with the a Ori observations (red dots with error bars) by iPerrin et al 
Top panel: intensity map in the K222-FLUOR filter of the snapshot that best matches the interferometric data. The range is 
[0;10 5 ] erg cm" 2 s" 1 A" 1 . The stellar parameters of this snapshot are: L = 93 480 L , 7? = 833 R , T eff = 3497 K and log g = -0.34. 
The simulation has been scaled to an apparent diameter of 43.6 mas at a distance of 179 pc. The white line indicates the position an- 
gle of the projected baselines. Other panels: Synthetic visibilities from the simulation, compared with observations. Panel B) covers 
the whole observational range, and panels C)-E) are close-ups of the clusters of observations. The thick solid line corresponds to 
the best match visibility curve(reduced^ 2 =0.21). The thin solid lines show the minimu m and maxim um extent of variations of the 
visibilities. The dot-dashed, and the dashed lines are the LD, and the UD models used by lPerrin et al.l (reduced^ 2 =22. 3, and 19.9). 
Note the logarithmic visibility scale. 
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